#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "mpi.h"
#include <memory.h>
#include "ZQ_CPP_LIB/time_cost.hpp"

#define PRINTR(format,...) \
if (rank==0) printf(format,##__VA_ARGS__)

#define PRINTI(format,...) \
printf("Proc#%d:\t" format,rank,##__VA_ARGS__)

#define PRINTRD(format,...) \
MPI_Barrier(MPI_COMM_WORLD); PRINTR(format,##__VA_ARGS__)

#define PRINTID(format,...) \
MPI_Barrier(MPI_COMM_WORLD); PRINTI(format,##__VA_ARGS__)


#define  MAX_PROCESSOR_NUM  32   /*set the max number of the Processor*/
#define  MAX_ARRAY_SIZE     1024*8   /*set the max size of the array*/
double a[MAX_ARRAY_SIZE][MAX_ARRAY_SIZE], g[MAX_ARRAY_SIZE][MAX_ARRAY_SIZE];


int main(int argc, char *argv[])
{
    int n;
    double transTime = 0,tempCurrentTime, beginTime;

    MPI_Status status;
    int rank, size;
    FILE *fin;
    int i, j, k;


    MPI_Init(&argc, &argv);
    MPI_Comm_rank(MPI_COMM_WORLD,&rank);
    MPI_Comm_size(MPI_COMM_WORLD,&size);
    
    n = atoi(argv[1]);


    auto t1 = zq_cpp_lib::get_timestamp();



    MPI_Bcast(&n, 1, MPI_INT, 0, MPI_COMM_WORLD);

    for(k = 0; k < n; k ++)
    {
        /* gathering the result ,and then broacasting to each processor */
        MPI_Bcast(a[k], (n-k)*MAX_ARRAY_SIZE, MPI_DOUBLE, 0, MPI_COMM_WORLD);
        PRINTI("line 85:\tk=%d\n",k);
                                                                                                            
        for(i = k+rank; i < n; i += size)
        {
            for(j = 0; j < k; j ++)
            {
                g[i][j] = a[i][j];
            }
            if (i == k)
            {
                for(j = k; j < n; j ++) g[i][j] = a[i][j]/sqrt(a[k][k]);
            }
            else
            {
                g[i][k] = a[i][k]/sqrt(a[k][k]);
                for(j = k+1; j < n; j ++) g[i][j] = a[i][j] - a[i][k]*a[k][j]/a[k][k];
            }
        }
        PRINTI("line102:\tk=%d\tReady to call MPI_Send\n",k);

        /* use the Cholersky Algorithm */
        if (rank!=0)
        for(i = k +rank; i < n; i ++)
        {
            
            PRINTI("k=%d\tSend:[%d] -> proc#%d\n",k,i,0);
            MPI_Send(g[i], n, MPI_DOUBLE, 0, k*1000+i, MPI_COMM_WORLD);
        }

        PRINTI("line112:\tk=%d\tsize=%d\tReady to call MPI_Recv\n",k,size);
        if(rank == 0)
        {
            // for (i = k; i < n; i++)
            //     memcpy(a[i],g[i],sizeof(double)*n);

            for(j = 1; j < size; j ++)
            {
                for(i = k + j; i < n; i += 1)
                {
                    PRINTI("k=%d\tRecv:[%d] <- proc#%d\n",k,i,j);
                    MPI_Recv(a[i], n, MPI_DOUBLE, j, k*1000+i, MPI_COMM_WORLD, &status);
                }
            }
        }
        PRINTI("line123:\tk=%d\tready for next iteration...\n",k);
    }
    PRINTI("OVER, to print results.\n");
    auto t2 = zq_cpp_lib::get_timestamp();
    PRINTR("%lf\n", zq_cpp_lib::get_cost_time_by_us(t1,t2));

    if (rank == 0)
    {
        puts("After Cholersky Discomposion");
        puts("Output Matrix G");
        for(i =0; i < n; i ++)
        {
            for(j = 0; j < i; j ++) printf("           ");
            for(j = i; j < n; j ++) printf("%9.5f  ", a[i][j]);
            printf("\n");
        } /* output the result */
    }

    MPI_Finalize();/* end of the program */
    return 0;
}

